function J = cost_function(X, U, e, data, P_load_future)
% =========================================================================
% cost_function.m - (修改版) NMPC的自定义成本函数
%
% 核心修改:
% 1. 更新了从状态矩阵X中提取 P_fc, SOC, alpha_decay 的索引，以匹配
%    新的状态向量定义 [SOC; T_stack; alpha_decay; P_fc_actual]。
% =========================================================================
% 输入:
%   X:      未来状态轨迹 [x(k+1|k); x(k+2|k); ...; x(k+p|k)]
%           维度: PredictionHorizon+1 x nx
%   U:      未来控制序列 [u(k|k);   u(k+1|k); ...; u(k+p-1|k)]
%           维度: PredictionHorizon x nu
%   P_load_future: 预测的未来扰动序列

    % 加载参数
    persistent CostParams SimParams
    if isempty(CostParams)
        CostParams = evalin('base', 'CostParams');
        SimParams  = evalin('base', 'SimParams');
    end
    
    dt = SimParams.Ts;
    J = 0; % 初始化总成本
    
    % --- 遍历整个预测时域，累加每一步的成本 ---
    for k = 1:data.PredictionHorizon
        % 提取当前预测步(k)和下一步(k+1)的状态
        state_k   = X(k, :);
        state_kp1 = X(k+1, :);
        
        % 根据新的状态向量定义更新索引
        SOC_k         = state_k(1);
        alpha_decay_k = state_k(3);
        
        SOC_kp1         = state_kp1(1);
        alpha_decay_kp1 = state_kp1(3);
        P_fc_kp1        = state_kp1(4); % P_fc 是第4个状态
        
        % =================================================================
        % 1. 氢气消耗成本 C_H2
        % =================================================================
        mass_h2_consumed = CostParams.k_H2 * P_fc_kp1 * dt;
        cost_h2 = CostParams.p_H2 * mass_h2_consumed;

        % =================================================================
        % 2. 锂电池寿命损耗成本 C_bat_loss
        % =================================================================
        delta_soc = abs(SOC_kp1 - SOC_k);
        cost_per_soc_change = CostParams.C_bat_total / CostParams.N_cycle;
        cost_bat = cost_per_soc_change * delta_soc;

        % =================================================================
        % 3. 燃料电池老化成本 C_fc_age
        % =================================================================
        delta_alpha_decay = alpha_decay_kp1 - alpha_decay_k;
        cost_per_alpha_change = CostParams.C_fc_total / CostParams.alpha_max;
        cost_fc = cost_per_alpha_change * delta_alpha_decay;
        
        % --- 累加总成本 ---
        J = J + cost_h2 + cost_bat + cost_fc;
    end
end
